Let’s import important libraries

library(dplyr)
library(sf)
library(leaflet)
library(ggplot2)
library(tibble)
library(rgdal)
library(proj4)
library(remotes)
remotes::install_github('ropensci/osmdata')
library(osmdata)
library(osmextract)
# proj_db <- system.file("proj/proj.db", package = "sf")

Import vector data from Open Street Map

# from : https://cran.r-project.org/web/packages/osmdata/vignettes/osmdata.html


bb <- getbb('rotterdam, NL', format_out = 'sf_polygon')
x <- opq(bbox = bb) %>% 
   add_osm_feature(key = 'building') %>%
    osmdata_sf ()

str(x$osm_polygons)
## Classes 'sf' and 'data.frame':   253568 obs. of  274 variables:
##  $ osm_id                               : chr  "23270254" "26545811" "34294334" "34294530" ...
##  $ name                                 : chr  NA "Montevideo" "Lotus Serre" "Ingang" ...
##  $ X3dr.height1                         : chr  NA NA NA NA ...
##  $ X3dr.length1                         : chr  NA NA NA NA ...
##  $ X3dr.type                            : chr  NA NA NA NA ...
##  $ X3dshapes.note                       : chr  NA NA NA NA ...
##  $ access                               : chr  NA NA "customers" NA ...
##  $ addr.city                            : chr  NA NA NA NA ...
##  $ addr.country                         : chr  NA NA NA NA ...
##  $ addr.full                            : chr  NA NA NA NA ...
##  $ addr.housename                       : chr  NA NA NA NA ...
##  $ addr.housenumber                     : chr  NA NA NA NA ...
##  $ addr.postcode                        : chr  NA NA NA NA ...
##  $ addr.street                          : chr  NA NA NA NA ...
##  $ addr.unit                            : chr  NA NA NA NA ...
##  $ aeroway                              : chr  NA NA NA NA ...
##  $ air_conditioning                     : chr  NA NA NA NA ...
##  $ alt_name                             : chr  NA NA NA NA ...
##  $ amenity                              : chr  NA NA "restaurant" NA ...
##  $ architect                            : chr  NA NA NA NA ...
##  $ area                                 : chr  NA NA NA NA ...
##  $ artist_name                          : chr  NA NA NA NA ...
##  $ artwork_type                         : chr  NA NA NA NA ...
##  $ atm                                  : chr  NA NA NA NA ...
##  $ bar                                  : chr  NA NA NA NA ...
##  $ bicycle                              : chr  NA NA NA NA ...
##  $ bicycle_parking                      : chr  NA NA NA NA ...
##  $ brand                                : chr  NA NA NA NA ...
##  $ brand.wikidata                       : chr  NA NA NA NA ...
##  $ brand.wikipedia                      : chr  NA NA NA NA ...
##  $ bridge.support                       : chr  NA NA NA NA ...
##  $ building                             : chr  "transportation" "apartments" "yes" "yes" ...
##  $ building.architecture                : chr  NA NA NA NA ...
##  $ building.colour                      : chr  NA NA NA NA ...
##  $ building.flats                       : chr  NA NA NA NA ...
##  $ building.levels                      : chr  "2" "43" NA NA ...
##  $ building.levels.underground          : chr  NA NA NA NA ...
##  $ building.material                    : chr  NA NA NA NA ...
##  $ building.min_level                   : chr  NA NA NA NA ...
##  $ building.part                        : chr  NA NA NA NA ...
##  $ building.use                         : chr  NA NA NA NA ...
##  $ capacity                             : chr  NA NA NA NA ...
##  $ changing_table                       : chr  NA NA NA NA ...
##  $ check_date                           : chr  NA NA NA NA ...
##  $ check_date.opening_hours             : chr  NA NA NA NA ...
##  $ cinema.3D                            : chr  NA NA NA NA ...
##  $ climbing.length.max                  : chr  NA NA NA NA ...
##  $ climbing.sport                       : chr  NA NA NA NA ...
##  $ clothes                              : chr  NA NA NA NA ...
##  $ club                                 : chr  NA NA NA NA ...
##  $ community_centre                     : chr  NA NA NA NA ...
##  $ company                              : chr  NA NA NA NA ...
##  $ construction                         : chr  NA NA NA NA ...
##  $ contact.website                      : chr  NA NA NA NA ...
##  $ content                              : chr  NA NA NA NA ...
##  $ covered                              : chr  NA NA NA NA ...
##  $ craft                                : chr  NA NA NA NA ...
##  $ cuisine                              : chr  NA NA NA NA ...
##  $ delivery                             : chr  NA NA NA NA ...
##  $ demolished                           : chr  NA NA NA NA ...
##  $ denomination                         : chr  NA NA NA NA ...
##  $ denomination.nl                      : chr  NA NA NA NA ...
##  $ description                          : chr  NA NA NA NA ...
##  $ description.de                       : chr  NA NA NA NA ...
##  $ designation                          : chr  NA NA NA NA ...
##  $ dhm_id                               : chr  NA NA NA NA ...
##  $ disused                              : chr  NA NA NA NA ...
##  $ dock                                 : chr  NA NA NA NA ...
##  $ drinking_water                       : chr  NA NA NA NA ...
##  $ email                                : chr  NA NA NA NA ...
##  $ emergency                            : chr  NA NA NA NA ...
##  $ facebook                             : chr  NA NA NA NA ...
##  $ fax                                  : chr  NA NA NA NA ...
##  $ fee                                  : chr  NA NA NA NA ...
##  $ ferry                                : chr  NA NA NA NA ...
##  $ fixme                                : chr  NA NA NA NA ...
##  $ floating                             : chr  NA NA NA NA ...
##  $ food                                 : chr  NA NA NA NA ...
##  $ foot                                 : chr  NA NA NA NA ...
##  $ frequency                            : chr  NA NA NA NA ...
##  $ fuel.GTL_diesel                      : chr  NA NA NA NA ...
##  $ fuel.diesel                          : chr  NA NA NA NA ...
##  $ fuel.lpg                             : chr  NA NA NA NA ...
##  $ fuel.octane_95                       : chr  NA NA NA NA ...
##  $ fuel.octane_98                       : chr  NA NA NA NA ...
##  $ generator.source                     : chr  NA NA NA NA ...
##  $ government                           : chr  NA NA NA NA ...
##  $ healthcare                           : chr  NA NA NA NA ...
##  $ healthcare.speciality                : chr  NA NA NA NA ...
##  $ height                               : chr  NA "152.3" NA NA ...
##  $ heritage                             : chr  NA NA NA NA ...
##  $ heritage.operator                    : chr  NA NA NA NA ...
##  $ historic                             : chr  NA NA NA NA ...
##  $ image                                : chr  NA "https://upload.wikimedia.org/wikipedia/commons/f/f2/Montevideo%28Rotterdam%29.jpg" NA NA ...
##  $ image.0                              : chr  NA NA NA NA ...
##  $ image.1                              : chr  NA NA NA NA ...
##  $ image.2                              : chr  NA NA NA NA ...
##  $ image.3                              : chr  NA NA NA NA ...
##  $ indoor                               : chr  NA NA NA NA ...
##   [list output truncated]
##  - attr(*, "sf_column")= chr "geometry"
##  - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA NA NA NA NA NA NA NA NA ...
##   ..- attr(*, "names")= chr [1:273] "osm_id" "name" "X3dr.height1" "X3dr.length1" ...

Map building age focusing on post-war buildings

buildings <- x$osm_polygons %>% st_transform(.,crs=28992)

buildings$build_date <- ifelse(buildings$start_date <1945, 1945,buildings$start_date)
 ggplot(data = buildings) +
   geom_sf(aes(fill = as.numeric(build_date), colour=as.numeric(build_date)))  +
   scale_fill_viridis_c(option = "viridis")+
   scale_colour_viridis_c(option = "viridis")

Let’s focus on really old building and imagine we’re in charge of conservation. We want to know how much of the city would be affected by a non-construction zone of 500m around pre-1800 buildings.

Selecting old buildings

old <- 1800 # year prior to which you consider a building old
 distance <- 500 # in meters
 
old_buildings <- filter(buildings, start_date <= old)
 ggplot(data = old_buildings) + geom_sf()

Creating a buffer around a polygon

 buffer_old_buildings <- st_buffer(old_buildings, dist = distance)
 ggplot(data = buffer_old_buildings) + geom_sf()

Fusing overlapping buffers into a single polygon

 single_old_buffer <- st_union(buffer_old_buildings) %>% 
   st_cast(to = "POLYGON") %>% st_as_sf()  %>% 
   add_column("ID"=as.factor(1:nrow(.)))  %>% st_transform(.,crs=4326) 

Representing old buildings by their centroid

sf::sf_use_s2(FALSE)
centroids_old <- st_centroid(old_buildings) %>% st_transform(.,crs=4326)  
ggplot() + 
    geom_sf(data = single_old_buffer, aes(fill=ID)) +
    geom_sf(data = centroids_old)

Counting the number of centroid per buffer with st_intersection (/! ≠st_intersect!)

 centroids_buffers <- st_intersection(centroids_old,single_old_buffer) %>% 
   add_column(n=1)
 centroid_by_buffer <- centroids_buffers %>% 
   group_by(ID) %>%
   summarise(
   n = sum(n)
   )
 single_buffer <- single_old_buffer %>% add_column(n_buildings = centroid_by_buffer$n)

 ggplot() + 
   geom_sf(data = single_buffer, aes(fill=n_buildings)) +
   scale_fill_viridis_c(alpha = 0.8,
                        begin = 0.6,
                        end = 1,
                        direction = -1,
                        option = "B")

Final visualisation

ggplot() + 
   geom_sf(data = buildings) +
   geom_sf(data = single_buffer, aes(fill=n_buildings), colour = NA) +
   scale_fill_viridis_c(alpha = 0.6,
                        begin = 0.6,
                        end = 1,
                        direction = -1,
                        option = "B") 

 old <- 1939 # year prior to which you consider a building old
 distance <- 100 # in meters

Now let’s look at the same analysis but for pre-war buildings and 100m buffers

 old <- 1939 
 distance <- 100
 #select
 old_buildings <- filter(buildings, start_date <= old)
 #buffer
  buffer_old_buildings <- st_buffer(old_buildings, dist = distance)
  #union
 single_old_buffer <- st_union(buffer_old_buildings) %>% 
   st_cast(to = "POLYGON") %>% st_as_sf()  %>% 
   add_column("ID"=1:nrow(.))  %>% st_transform(.,crs=4326) 
 #centroids
 centroids_old <- st_centroid(old_buildings) %>% st_transform(.,crs=4326)  
  #intersection
  centroids_buffers <- st_intersection(centroids_old,single_old_buffer) %>% 
   add_column(n=1)
 centroid_by_buffer <- centroids_buffers %>% 
   group_by(ID) %>%
   summarise(
   n = sum(n)
   )
 single_buffer <- single_old_buffer %>% add_column(n_buildings = centroid_by_buffer$n)
  ggplot() + 
   geom_sf(data = buildings) +
   geom_sf(data = single_buffer, aes(fill=n_buildings), colour = NA) +
   scale_fill_viridis_c(alpha = 0.6,
                        begin = 0.6,
                        end = 1,
                        direction = -1,
                        option = "B") 

Problem: there are many pre-war buildings and the buffers are large so the number of old buildings is not very meaningful. Let’s compute the density of old buildings per buffer zone.

Visualisation

single_buffer$area <- st_area(single_buffer, ) %>% units::set_units(., km^2)
 single_buffer$old_buildings_per_km2 <- as.numeric(single_buffer$n_buildings / single_buffer$area)

 ggplot() + 
   geom_sf(data = buildings) +
   geom_sf(data = single_buffer, aes(fill=old_buildings_per_km2), colour = NA) +
   scale_fill_viridis_c(alpha = 0.6,
                        begin = 0.6,
                        end = 1,
                        direction = -1,
                        option = "B") 

EXERCISE!

Get OSM limits of Rotterdam

 rdam <- getbb(
  place_name = "Rotterdam",
  format_out = 'sf_polygon',
  silent = FALSE
) 
## [1] "https://nominatim.openstreetmap.org/search?format=json&q=Rotterdam&featuretype=settlement&polygon_text=1&limit=10"

Import and reproject CBS 500mx500m grid cells for the entire NL

# downloaded from: https://www.cbs.nl/-/media/cbs/dossiers/nederland-regionaal/vierkanten/500/2022-cbs_vk500_2019_vol.zip 
cells_NL <- st_read("2019-cbs_vk500_2015_v2/CBS_VK500_2015_v2.shp") %>%
  mutate_if(is.numeric, funs(ifelse(.==-99997, NA, .)))
## Reading layer `CBS_VK500_2015_v2' from data source 
##   `/Users/ccottineau/Documents/Teaching/DataCarpentryWorkshop/GeoSpatial/2019-cbs_vk500_2015_v2/CBS_VK500_2015_v2.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 151108 features and 131 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 13000 ymin: 306500 xmax: 278500 ymax: 619500
## Projected CRS: Amersfoort / RD New
cells_NL <- st_transform(cells_NL,crs=4326) 

INTERSECT (≠ st_intersection) to keep only cells for rotterdam

intersect_cells_rdam <- st_intersects(cells_NL, rdam, sparse=F)
intersect_cells_rdam_sf <- cbind(intersect_cells_rdam,cells_NL)
cells_rdam <- intersect_cells_rdam_sf %>% 
  filter(X1 == T ) 

plot(cells_rdam)